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dispersive plasmonic nanostructures and metamaterials. We specifically apply it 
to dielectric refractive index, and metallic plasma frequency modulation in metal- 
dielectric nanostructures. As a numerical demonstration, we verify the theory's accu- 
racy against direct calculations, for a system of plasmonic rods in air where the metal 



is defined by a two-pole fit of silver's dielectric function. We also discuss new optical 
behavior related to plasma frequency modulation in such systems. Our approach 



provides new physical insight for the design of plasmonic devices for biochemical 
sensing and optical modulation, and future active metamaterial applications. 
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To design active optical devices such as sensors, switches and modulators, one needs 
to calculate how a small change in refractive index affects the device response function. 
For active devices based on dielectric structures^"- described by a frequency-independent 
dielectric constant distribution e(v), the effect of an index change can be understood in 
terms of a frequency shift oj\ of the eigenmodes of the system, which is given by first-order 
perturbation theory as:-£ 

q; /^A g (r)|E(r)| 2 
Ul "~ 2 / rfrc(r)|E(r)| 2 1 ' 

The numerator in Eq. (1) only has contributions from the perturbed regions as described 

by Ae(r). The shift in the eigenfrequency thus depends on the overlap of the modal electric 

field energy with the perturbed region. 

In recent years, there has been substantial interest in using plasmonics for active devices. 

While surface-plasmon sensors are already prominent in biochemical sensing applications,- 1 ^ 
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remarkable improvements in device performance have been achieved using plasmonic nanos- 
tructures and metamaterials.-~— Active plasmonic devices have also been implemented for 
modulation and switching.— ~— Given these developments, it is essential to develop a general 
photonic perturbation theory for plasmonic nanophotonic structures to firmly ground the 
analysis and characterization of next-generation devices. 

Eq. (1), however, is not applicable for plasmonic systems. For example, e(r) can be 
negative in a metal system, and hence directly applying Eq. (1) could lead to a prediction of 
infinite sensitivity, which is unphysical. Moreover, in plasmonic systems, recent experiments 
have successfully varied the plasma frequency of the metal as well,— - — introducing a new 
degree of freedom that requires formal theoretical treatment. 

In this paper, we build upon a formulation recently developed to model the photonic 
band structure of dispersive material systems, 20 to construct a perturbation theory that 
can predict modal frequency shifts due to changes in the dielectric constants of dispersive 
systems involving both metals and dielectrics. We also show that a similar perturbation 
theory can be developed when the plasma frequency of the dispersive metal is modulated, 
and highlight new physical behavior related to such modulation. Our theory thus allows us 
to treat variations in key parameters of both metal and dielectric components of plasmonic 
nanostructures. 

The paper is organized as follows. In Section II we review the photonic band theory 
previously developed to describe dispersive metamaterials and plasmonic nanostructures. 
We develop a general perturbation theory in Section III to describe the effect of small 
parameter changes on the optical modes of dispersive nanostructures. We consider the 
specific cases of refractive-index variation in Section IV, and plasma frequency variation in 
Section V, and demonstrate the accuracy of this perturbation theory numerically for both 
cases by comparing to direct calculations. Finally, we conclude in Section VI. 

II. REVIEW OF DISPERSIVE METAMATERIAL BAND THEORY 

As the basis for the developments of this paper, we first briefly review the theory recently 
developed to model the photonic modes of dispersive material systems.— This theory is 
applicable for a nanophotonic structure containing a dispersive material with a dielectric 
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function 

e(w) = £00 + £00 V — Wp 2 ' n , . r ■ (2) 

To describe such a dispersive material, for the n-th pole in the dielectric function, one 
introduces a polarization field P n and a polarization velocity field V n , satisfying 

BP 

^ = V„ (3) 
<9V 

—gjT = w p,n £ ooE - W n P n - T n V n (4) 

These auxiliary fields 2 ^ 2 ^ are then coupled to Maxwell's equations through: 

f^VxE (5) 

at no 



dE 

~dt ~ e 



1 N 

-(VxH-^V„) (6) 



n=l 



For steady state, with fields varying as exp(iojt), Eqs. (3)- (6) become 

iu~H = -— V x E (7) 
Mo 



<wE = -(VxH-^V n ) (8) 

n=l 



iwP n = V n (9) 
iwV n = Wp in £ooE - wJ n P„ - T n V„ (10) 

and thus define an eigenvalue problem for to. Eqs. (jTjl- tfTOj) also define a total energy density 

1 N 1 

W = -(5oo|E| 2 + /i |H| 2 ) + TT-^KnlPnl 2 + |V n | 2 ) (11) 

n=l 

The spatial integral of Eq. ( II ip . which represents the total energy of the system, is conserved 
when T n = for all poles. Further, defining x = (H, E, Pi, Vi, ■ • • ,Ptv,Vtv) T , which 
represents a multi-component vector field varying over the whole space, we can write Eqs. 

(ED-dloD as 

wAx = Bx, (12) 
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where A = diag (fi , £«,,••-, w 2 ,iv/ w p,iv £ oo, and 
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(13) 
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Finally, defining y = v Ax, we can re-write this as an eigenvalue equation for u: 

W y=(v / A)" 1 B(v / A)" 1 y. (14) 

For the lossless case, T n = 0, Eq. ( IT41 becomes a Hermitian eigenvalue equation which 
results in an orthogonality condition: 

N 



[ dr \( £oo ei • E b + /i H; • H b ) + £ ^^ Wo 2 n p; in • p b , n + v; n • v b , r 



Sab 

(15) 



III. PERTURBATION THEORY 



In this section, we develop a general perturbation theory based on the generalized eigen- 
value equation for the photonic bands of dispersive nanostructures, Eq. ( Tl2l) . We start from 
the unperturbed system 

a; A xo = B x (16) 

In the presence of a perturbation, the system matrices become A = Ao + Ai, and B = 
B + Bi, and result, we have 

(u Q + Wi)(Aq + Ax)(xo + Xi) = (B + Bi)(xo + x x ) (17) 

Using Eq. (1161) and keeping only first-order terms in Eq. (j!7p . we have 

w A xi + cjiA x + w Aix = B xi + Bix (18) 

Eq. ffTBT) . in its most general form, describes a lossy system and cannot be written as a 
Hermitian eigenvalue problem. Thus, to calculate uj\ we also need to determine the left 
eigenvector z that satisfies 

a; z Ao = z B (19) 



5 



Multiplying z through Eq. (1TB]) and solving for ui, we find: 

z Bix - w z Aix 

= (20) 

z A x 

Eq. (120]) is the main result of this paper. In the following sections, we apply this equation 
to two relevant examples of perturbations in plasmonic systems. 



IV. APPLICATION I: DIELECTRIC REFRACTIVE-INDEX MODULATION 

In this section we consider the specific case of a small change (Ae(r)) in the dielectric 
constant of a dielectric region, in a nanostructure consisting of both metal and dielectric 
regions. The metal region is assumed to be unperturbed. In this case, the perturbation 
takes the form 

A 1 = diag(0,Ae(r),-.. ,0,0), (21) 

while Bi = 0. We now determine the change in modal frequency for the cases when the 
metal in the metal-dielectric nanostructure is lossless and lossy. 



A. Lossless case, T n = 

For the lossless case Eq. ( 112]) both A and B are Hermitian. In this case, from Eq. ()19p . 
we have z = xj. Thus, Eq. (120]) reduces to 

xjAiXp J^rA g (r)|E(r)| 2 

wi = -Uq-t- = -w f — . (22) 

x A x / dr Wo 

For the lossless dispersive system, we thus obtain a result that has the same form as Eq. (1) 
that is now appropriate for a system with dispersion, provided that we consider the total 
energy density in the system including contributions from the auxiliary mechanical fields. 
The expression for total energy, Eq. ( TTTT) . has contributions from multiple Lorentz poles 
and is a multi-pole extension of the energy density expression previously derived by taking 
electric polarization into account explicitly.— 1 ^ For the lossless case this reduces to the usual 
expression for energy density in metals 23 



W = \ 



d(0JEE c 

du 



|E| 2 + ^|H| 2 . (23) 
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B. Lossy case, T n ^ 

For the lossy case, the matrix B in Eq. (fTBT) is no longer Hermitian and z ^ xj. Thus, 
Eq. ( 1201) in this case reduces to: 

ZoAiXq 

= — t (24) 

z A x 



While no explicit expression analogous to Eq. (1) can be written for the lossy case, Eq. ((24 
still allows one to calculate the frequency shift due to a dielectric refractive index change in 
the presence of a lossy metal; an important ability in realistic plasmonic sensing schemes. 

Eq. f)24p represents the technically correct way to do perturbation theory, where one needs 
to determine both the left and right eigenvectors of the general eigenvalue problem. Moreover 
the denominator in Eq. ([53]) cannot be interpreted as an energy integral. Empirically, on 
the other hand, we will show numerically that in fact zq ps xj, even for metals with realistic 
loss parameters, and thus the denominator of Eq. (124]) . 

z A x « J drW (25) 

where Wo is the energy density of the mode for the lossy system as defined in Eq. (llip . 
which includes contributions from the mechanical auxiliary fields. We note that in a lossy 
system, when multiple poles are involved, there is no simple relation such as Eq. (123]) that 
can be used to describe the total energy. Instead, the definition of Eq. ( ITT]) , which explicitly 
takes into account contributions from the auxiliary mechanical fields, must be used. 

C. Numerical Example 

We now numerically verify the accuracy of the perturbation theory results presented. 
Motivated by a recent experiment,- we consider as our model system a two-dimensional 
periodic array of square plasmonic rods in air that are uniform along the third z direction, 
and numerically solve its eigenmodes. The system has periodicity a = 130nm and the rod 
has a side length of s = 0.45a = 58.5nm. 

1. Metal Dielectric Function Fit 

The plasmonic metal's dielectric function is defined by fitting silver's tabulated 
permittivity 25 with a Lorentz-Drude model,— consisting of a Lorentz pole and a Drude 
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FIG. 1: (a) Comparing tabulated data for the real part of silver's permittivity 25 against the Drude 
model and Lorentz-Drude 2 - fit used in this paper, at optical wavelengths, (b) The imaginary 
and real parts of the Lorentz-Drude fit 2 ^ of silver's permittivity over normalized frequencies for 
a = 130nm. 

pole [Eq. (2)]. In Fig. 1(a) we see that using the two poles is more accurate than using 
the Drude model alone for modeling silver's dielectric function at optical frequencies. The 
Lorentz pole in this fit is denned by w ,i = 0.5526, co p ^ = 0.8196, and T 1 = 0.1195, and the 
Drude pole by u 0:2 = 0, uj p>2 = 0.9615, and T 2 = 0.0022. All frequencies in the numerical 
examples are normalized to 27rc/a. A plot of the real and imaginary parts of this fit over the 
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FIG. 2: Computed band structures for a square lattice of 2D plasmonic rods (s = 0.45a) in air 
for the (a) TM and (b) TE polarization between (k x = 0, k y = 0) and (k x = Tr/a,k y = 0) points. 
The metal is describe by a Lorentz-Drude fit of silver with the presence of the Lorentz pole at u)q \ 
highlighted. The cutoff frequency for TM modes is identified at u c , and frequency regions where 
Re[e] > are highlighted in yellow. 

relevant normalized frequency range is presented in Fig. 1(b). In the formalism of Section 
II, the Drude pole has zero resonant frequency, and hence only requires the V-field as its 
auxiliary field. Thus, in our system, we describe the effects of dispersion in terms of three 
auxiliary fields: Pi and Vi for the Lorentz pole, and V 2 for the Drude pole. 



2. Theory Verification 

We analyze the TM and TE modes (with their electric and magnetic fields respectively 
polarized along the ^-dimension) of this system in Fig. 2. These band structures are 
calculated by implementing a finite-difference spatial discretization of the fields with a Yee 
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FIG. 3: Comparing the perturbation theory prediction ('Perturbed') and direct solution ('Direct') 
of uj\ in the region identified in Fig. 2(b), for dielectric Ae(r) = 0.02. The (a) lossless metal 
(r n = 0) and (b) lossy metal cases shows excellent agreement for both the flat surface mode and 
non-surface mode. 

grid, and then solving the eigenvalue equation using the Arnoldi method.— In our finite 
difference implementation, we have truncated the finite-difference grid appropriately at the 
metal-air interface to ensure that boundary conditions for the tangential field components 
are satisfied at these interfaces.— 1 ^ For the TM case, we find that a stop-band exists below 
a cutoff frequency u c = 0.3067. In the TE case we note the presence of dispersion-less 
flat bands below u> ~ 0.4 that correspond to surface plasmon modes, and modes that have 
substantial group velocity which correspond to non-surface modes. For both cases, we also 
observe the presence of a dense cluster of low group-velocity modes right below cuo,i— in the 
region where Re[e] > 0, highlighted yellow in Fig. 1(b) and Fig. 2. 

To verify the accuracy of our perturbation theory approach we highlight a region of the 
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FIG. 4: Reflectivity spectrum of a finite 2D square lattice of square plasmonic rods (s = 0.45a) in 
air, consisting of 50 layers. The results are obtained with a full-field FDFD simulation (shown in 
the inset schematic) and show a reflectivity dip corresponding to the lowest-frequency propagating 
mode in the system (the mode at uj c in Fig. 2(a)). Altering the air region by Ae = le — 3 to 
simulate a perturbation, a shift in the dip of Ato = — 9.9e — 5 is observed, matching the theoretical 
prediction of Eq. ([22]) . uj\ = — 9.37e — 5, well. 

TE band structure featuring both a surface and non-surface mode in Fig. 3. We alter the 
dielectric constant of the dielectric region by Ae = 0.02, and calculate the resulting shift in 
eigenfrequency, using perturbation theory for both lossy and lossless cases. For the lossless 
system, we set T n = 0, and use Eq. ( I2"2"j) . For the lossy system, we can use the exact 
perturbation theory result of Eq. (124"]) . as well as the approximation in terms of energy 
density in Eq. fl25|) . The results from these two forms of perturbation theory are nearly 
identical to each other. These results, from both forms of the perturbation theory, are then 
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compared to the band structure obtained by directly solving Eq. ( JT2I) for the perturbed 
system. The results from the perturbation theory show excellent agreement with results 
from the direct calculation for both the lossless and lossy systems. 

3. Practical Application 

As an illustration of the practical significance of this result, we analyze a potential 
refractive- index sensing scheme by calculating the reflection/transmission spectrum of 50- 
layers of the plasmonic rod system considered above, using a full-field 2D finite-difference 
frequency- domain (FDFD) simulation. 30 The plasmonic metal, corresponding as before to 
the two-pole fit of silver's dielectric function, is assumed to be lossless for the purposes of 
this illustration. For the TM case, below the cutoff frequency u c [Fig. 1(a)], no propagating 
modes are supported. Thus, the structure is strongly reflecting. The lowest-frequency dip 
in the reflection spectrum [Fig. 4] corresponds to the first mode supported by the system, 
identified previously in Fig. 2(a) at u c . We can then examine the shift in this reflectivity 
dip when the dielectric surrounding the rods is changed by Ae (due to, for example, the 
introduction of a biochemical analyte). We observe a shift of the dip by Au = — 9.9e — 5 
for Ae = le — 3. Using Eq. (|22|) we calculate the shift of the lowest-frequency mode at 
k = 0, for the same Ae, to be uji = — 9.37e — 5, which matches well with the shift observed 
in the full-field simulation. Thus, the perturbation theory with respect to the eigenmodes of 
the system can be used to predict shifts in features of the same system's transmission and 
reflection spectra. 

V. APPLICATION II: METALLIC PLASMA FREQUENCY MODULATION 

In this section we demonstrate that we can also treat small changes in the plasma fre- 
quency of a metal present in a plasmonic nanostructure. This is of interest given a recent 
experimental work with ITO^ where an applied electric field induced a change in carrier 
concentration thereby changing its plasma frequency and behavior at optical frequencies. In 
other recent work, the infrared plasmonic resonance of a metallic nanostructure (split-ring 
resonators) immersed in an electrolyte was affected by applying a voltage that altered the 
structure's average electron bulk density— 
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To derive an expression for the modal frequency change in such situations, without loss 
of generality, we will assume that only the N-th pole experiences a plasma frequency mod- 
ulation. Then, we note that to first order the effect of Ao; P) jv affects our matrix elements by 
ojp ]iV — > ojp^ N + 2oj P! nAu> p ,n- With this established we now separately consider the lossless 
and lossy cases. 



A. Lossless case, T r 



For the lossless case the following perturbing matrices, Ai and Bi are applicable: 



Ai = diag 0, 0, 



0,N 



/2Aw. 



p,N 



U p,N £ oo V U P,N 



U p,N £ oo 



2 Aw 



p,N 



Up.N 



(26) 



••• 

••• 








( 2Aoj 



i'~ 2 

W N Ecx> 



2Ai 



>p,N 



>p,N 



'p.N 6x > \ ^P.N 







(27) 



/ 



Since the eigenvalue equation Eq. (TTzjl is Hermitian for the lossless case, Eq. ( 1201) becomes: 

xJBiXq - WqxJAiXq 



XqA Xo 



(28) 



which reduces to the following expression in terms of the constituent fields (as determined 
by solving Eq. f fT2|) ): 



OJQ 



J dr W 



dr 



2Aup, N (r) uJq-uj, 



0,N 



^p.N 



UJ, 



|Piv(r) 



(29) 



p,N 



The frequency shift uj\ for a given mode is thus directly proportional to the strength 
of the modal mechanical field in the iV-th pole that is undergoing modulation, and the 
proximity of the unperturbed modal frequency uq to uo t N- Of particular interest is that, 
assuming Au Pt N > 0, u± is positive or negative depending on whether the unperturbed 
modal frequency is above or below the Lorentz pole's resonance frequency w^tv- We verify 
this behavior numerically in Section C below. 
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B. Lossy case, T n ^ 



As noted in the previous section, for the lossy case, Eq. (TT2j) is no longer Hermitian and 
z 7^ xj. Thus, Eq. ( 1201) in this case reduces to: 

Z BiX - WqZoAiXo 



Ul z A x 
We note here that Bi is slightly altered in the lossy case: 

/o ••• 



Bi 



\ 



•• 
•• 







,2 



,2 

"0,N 



2A, 



IjV 



2Al 



>p,N 



( 2Aoj. 



'p,N 



~l 2 — I 



(30) 



(31) 



As in Sec. IV. B, we note that empirically at optical frequencies zo ~ xj, and the denom- 
inator of Eq. ~ / dr W . 



C. Numerical Example 

We consider the same example nanostructure of plasmonic rods in air used in the previous 
section, this time to numerically verify the accuracy of the perturbation theory for plasma 
frequency modulation in the plasmonic rod. Specifically, we focus again on the region of 
the TE band structure highlighted in Fig. 2(b) featuring both a surface and non-surface 
mode. As before, in addition to the original structure's TE bands, we directly calculate 
the bands of the perturbed structure, where the metal rod's Lorentz pole is modulated as 
Au Pt i = 0.050^1. We now use the perturbation theory calculations to predict the effect of 
Aw Pi i and compare it to the direct result. We note excellent agreement in Fig. 5 between 
the perturbation theory and direct approaches for both lossless and lossy cases. 

Moreover, we note that in both lossy and lossless cases, the surface plasmon mode expe- 
riences a substantially greater oj\ compared to the lower frequency non-surface mode. This 
is consistent with the observation from Raman and Fan22 that surface modes have stronger 
mechanical field intensities, and the form of Eqs. (129]) and fl30l . which predict greater U\ 
for modes with stronger mechanical fields. This behavior is of practical interest for active 
plasmonic device applications. 
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FIG. 5: Comparing the perturbation theory prediction ('Perturbed') and direct solution ('Direct') 
of u\ in the same zoomed-in region identified in Fig. 2(b), for Au Pi i = 0.05w Pi i. The (a) lossless 
(r n = 0) and (b) lossy cases shows excellent agreement for both the flat surface mode and non- 
surface mode. Note the substantially greater uj\ for the flat plasmon mode which has a stronger 
mechanical component to its eigenmode, and that u)\ < as predicted by Eq. (|29p . 

Finally, as observed in Eq. f )29|) . for a positive shift of the plasma frequency Aw Pjl > 0, 
we expect either a positive or negative frequency shift u\ depending on whether a mode's 
cu is greater or less than co>o,i- In Fig. 5 the modes presented are below wo,i and u\ < 0, as 
expected for Ao; Pi i > 0. In Fig. 6, we examine modes in the region above highlighted 
in Fig. 2(b), and note that uj x > as predicted by Eq. ( |29|) . This unusual property is 
highlighted here for the first time in the literature, and of potential utility as research in 
active plasmonics progresses. 
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FIG. 6: Comparing the perturbation theory prediction ('Perturbed') and direct solution ('Direct') 
of oj\ in the region above wo,i identified in Fig. 2(b), for Aw P) i = 0.05a/ Pj i. The (a) lossless (r n = 0) 
and (b) lossy cases show uj\ > (as opposed to lo\ < in Fig. 5), verifying the prediction of Eq. 

(22D. 

VI. CONCLUSION 

We have developed and verified a perturbation theory that provides considerable insight 
into the effect of small variations in both the dielectric and metallic components of plasmonic 
and dispersive metamaterial nanostructures. Using a previously developed photonic band 
formalism for general dispersive nanostructures described by an arbitrary number of Lorentz 
poles, we are able to accurately predict the effect of modulation in the dielectric refractive 
index and metallic plasma frequency of a metal-dielectric plasmonic system. We specifi- 
cally highlight new behaviors related to plasma frequency modulation that are of potential 
significance for future research in active plasmonics. 
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Given increasing interest in such dispersive systems for sensing and modulation appli- 
cations, and active devices more generally, the theory developed above will indeed be of 
considerable utility in future designs and analyses. 
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